
function [p q] = implicit_theta(dt,beta)
global p  q  nt x  L1 L;


I = eye(2*x);
% Periodic Boundary Conditions
sol(1:x,1) =p;
sol(x+1:2*x,1) =q;
%p_global(:,1)=p;

dt2= 1/dt;
for n=2:nt+1;
    r=L*sol*dt2;
    rhs = -[r;r];
        
    %rhs = (dt2*I - (1-beta)*L)*sol;
    %sol=(dt2*I+beta*L)\rhs;    
    K = L1\rhs;
    sol = sol+dt*0.5*(K(1:2*x,1)+K(2*x+1:4*x,1));
    p = sol(1:x,1);
    q= sol(x+1:2*x,1);
%p_global(:,n)=p;
    if rem(n,100)==0
        
        %refreshdata(h1,'caller') % Evaluate p in the function workspace
     %   drawnow
    end
    
end
display('Completed Successfully');